%Create 
FileNamesOuter=cell(1,6);
FileNamesOuter{1,1} = {'output\FirstWorkspace6AH20','output\SecondWorkspace6AM20','output\ThirdWorkspace6AU20','output\FourthWorkspace6AZ20','output\FifthWorkspace6AH21'};
FileNamesOuter{1,2} = {'output\FirstWorkspace6BH20','output\SecondWorkspace6BM20','output\ThirdWorkspace6BU20','output\FourthWorkspace6BZ20','output\FifthWorkspace6BH21'};
FileNamesOuter{1,3} = {'output\FirstWorkspace6CH20','output\SecondWorkspace6CM20','output\ThirdWorkspace6CU20','output\FourthWorkspace6CZ20','output\FifthWorkspace6CH21'};
FileNamesOuter{1,4} = {'output\FirstWorkspace6EH20','output\SecondWorkspace6EM20','output\ThirdWorkspace6EU20','output\FourthWorkspace6EZ20','output\FifthWorkspace6EH21'};
FileNamesOuter{1,5} = {'output\FirstWorkspace6JH20','output\SecondWorkspace6JM20','output\ThirdWorkspace6JU20','output\FourthWorkspace6JZ20','output\FifthWorkspace6JH21'};    
FileNamesOuter{1,6} = {'output\FirstWorkspace6SH20','output\SecondWorkspace6SM20','output\ThirdWorkspace6SU20','output\FourthWorkspace6SZ20','output\FifthWorkspace6SH21'};
%Create string with root letters for different markets, A=AUD/USD,
%B=GBP/USD, C=CAD/USD, E=EUR/USD, J=JPY/USD,S=CHF/USD
Strcatroots = {'A','B','C','E','J','S'};

%Change directories
%functionsFolder=cd('..');
%parentDirectory=cd('output\');
betaValues = (0.3:0.025:0.6)'; %13x1 vector of exponent values

nVec = [1];
outerloop = 0;
%outerloop cycles through different roots (contract markets)
%can change outerloop to 6 from 5 to run CHF/USD test
while(outerloop<5)
    tic
    outerloop = outerloop + 1;
    FileNames = FileNamesOuter{1,outerloop};
    n = nVec(1);

    %Initialize values for each storage variable
    
    %Beta Coefficients for different models and standard errors
    BetasMech = cell(11,1);
    BetasMech05d = cell(11,size(betaValues,1));
    BetasNMech = cell(11,1);
    StandardErrorsMech = cell(11,1);
    StandardErrorsMech05d = cell(11,size(betaValues,1));
    StandardErrorsNMech = cell(11,1);
    
    %Define general nonlinear model to estimate
    nlModel = @(Beta,X)(Beta(1)+Beta(2)*((abs(X).^Beta(3)).*sign(X))); 
    %Define square root price impact model
    nlModel2 = @(Beta,X)(Beta(1)+Beta(2)*((abs(X).^0.3).*sign(X))); 
    

    
    LambdasNMechL = cell(11,1);
    BetasNMechL = cell(11,1);
    SandwichesNMechL = cell(11,1);
    DaySizesNMechL = cell(11,1);
    R2sNMechL = cell(11,1);

    LambdasMechL = cell(11,1);
    BetasMechL = cell(11,1);
    SandwichesMechL = cell(11,1);
    DaySizesMechL = cell(11,1);
    R2sMechL = cell(11,1);
    VolumesMechL = cell(11,1);
    AdditionalValues = cell(11,1);

    DatesForLater = cell(11,1);

    %Four vectors to build out: price with mechanical depletion, price with
    %only non-mechanical, signed order size, date/time (Y1Mech Y1NMech X1 t2)
    Y1MechAll = [];
    Y1NMechAll = [];
    X1All = [];
    t2All = [];
    retsAll = [];

    filecounter = 0;
    %Cycle through each contract month per market: Mar '20, June '20, Sep
    %'20, Dec '20, Mar '21
    while(filecounter<5)
        filecounter = filecounter + 1;
        clear ESU1820180801;
        
        %Load table of data for contract month
        load(FileNames{1,filecounter});
        
        %Convert timezome from UTC to NYC time
        a2 = datetime(ESU1820180801{:,1},'TimeZone','UTC');
        a2.TimeZone = 'America/New_York';
        t2 = datenum(a2);
        
        %Clean up NaNs if any
        dataIsNaNFormat = double(ESU1820180801{:,2:7});
        dataIsNaNFormat(dataIsNaNFormat == 0)= NaN;
        
        %Run routine that analyzes dataset to create a list of aggressive
        %orders signed by direction, with size and order book dynamics
        %corresponding to order time
        [AggressiveOrders2, ~,Y2] = MicroStructureES([t2 dataIsNaNFormat]);
        
        %Add a column that consists of the "true price"
        AggressiveOrders2(:,13) = (AggressiveOrders2(:,8).*AggressiveOrders2(:,9) + AggressiveOrders2(:,10).*AggressiveOrders2(:,7))./(AggressiveOrders2(:,8) + AggressiveOrders2(:,10));
        
        %Filter out trades outside of liquid trading hours
        t2H = hour(AggressiveOrders2(:,1));
        t2M = minute(AggressiveOrders2(:,1));
        
        selector = ((t2H>=18) +(t2H<15)); %FX & US Rates
        %selector = ((t2H>=3) +(t2H<12)); %Europe Rates
        AggressiveOrders3 = AggressiveOrders2(selector>0.01,:);
        clear AggressiveOrders2;
        
        %Create X and Y vectors using careful elimination of intraday jumps
        t2H = hour(AggressiveOrders3(:,1));
        diffjumps = ((abs(t2H(2:end)-t2H(1:end-1))>2)+ (abs(t2H(2:end)-t2H(1:end-1))<22))>1.01;
        diffjumps = [diffjumps; 0]; %#ok<AGROW> %next order is following day
        
        %Determine number of days in the dataset
        selectorDays = diffjumps.*(1:size(diffjumps,1))';
        selectorDays = [0;selectorDays(selectorDays>0,:);size(AggressiveOrders3,1)];
        
        %Pre-allocate based on number of days for daily estimates
        betasMech = zeros(size(selectorDays,1)-1,3);
        betasMech05d = zeros(size(selectorDays,1)-1,2);
        betasNMech = zeros(size(selectorDays,1)-1,3);
        standardErrorsMech = zeros(size(selectorDays,1)-1,3);
        standardErrorsMech05d = zeros(size(selectorDays,1)-1,2);
        standardErrorsNMech = zeros(size(selectorDays,1)-1,3);
        lambdasNMech = zeros(size(selectorDays,1)-1,1);
        betasNMechL = zeros(size(selectorDays,1)-1,2);
        sandwichesNMech = zeros(size(selectorDays,1)-1,4);
        daysizesNMech = zeros(size(selectorDays,1)-1,1);
        r2sNMech = zeros(size(selectorDays,1)-1,1);

        lambdasMechL = zeros(size(selectorDays,1)-1,1);
        betasMechL = zeros(size(selectorDays,1)-1,2);
        sandwichesMechL = zeros(size(selectorDays,1)-1,4);
        daysizesMechL = zeros(size(selectorDays,1)-1,1);
        r2sMechL = zeros(size(selectorDays,1)-1,1);
        volumesMechL = zeros(size(selectorDays,1)-1,1);
        additionalValues = zeros(size(selectorDays,1)-1,6);
        
        datesForLater = zeros(size(selectorDays,1)-1,1);
        
        %cycle through beta values
        betacounter = 0;
        
        while(betacounter<size(betaValues,1))
            betacounter = betacounter+1;
            counter = 0;
            while(counter<size(selectorDays,1)-1)
                %Loop through days and run analysis
                counter = counter + 1;
                startval = selectorDays(counter,1)+1;
                endval = selectorDays(counter+1,1);
                datesForLater(counter,1) = AggressiveOrders3(startval+1,1);
                %n=1; %n controls the number of transactions ahead for which we
                %are calculating price impact
                %X1 is signed order size
                X1 = AggressiveOrders3(startval:endval-n,3).*AggressiveOrders3(startval:endval-n,4);
                %Y1Mech is the difference in true price from before one trade
                %to the subsequent trade
                Y1Mech = AggressiveOrders3(startval+n:endval,6)-AggressiveOrders3(startval:endval-n,6);
                %Y1NMech is the non-mechanical change in true price (price
                %immediately before execution minus the price immediately after execution of the previous aggressive order
                %Now square root estimate
                nlModel2 = @(Beta,X)(Beta(1)+Beta(2)*((abs(X).^(betaValues(betacounter,1))).*sign(X)));
                [betaMech05d,R105d,J05d,~,MSEMech05d,~] = nlinfit(X1,Y1Mech,nlModel2,[0;1]);
                CovMtrxMech05d = NLLSVarMtrxCalcs(J05d,R105d);
                CovMtrxMech05d = CovMtrxMech05d/size(R105d,1);
                betasMech05d(counter,:) = betaMech05d';
                standardErrorsMech05d(counter,:) = (diag(CovMtrxMech05d).^0.5)';
                betaMech05d'
            end
            BetasMech05d{filecounter,betacounter} = betasMech05d;
            StandardErrorsMech05d{filecounter,betacounter} = standardErrorsMech05d;
        end
        
        %Store values
        Y1MechAll = [Y1MechAll;AggressiveOrders3(1+n:end,6)-AggressiveOrders3(1:end-n,6)]; %#ok<AGROW>
        Y1NMechAll = [Y1NMechAll;AggressiveOrders3(1+n:end,6)-AggressiveOrders3(1:end-n,13)];%#ok<AGROW>
        X1All = [X1All; AggressiveOrders3(1:end-n,3).*AggressiveOrders3(1:end-n,4)];%#ok<AGROW>
        t2All = [t2All; AggressiveOrders3(1:end-n,1)]; %#ok<AGROW>
        retsAll = [retsAll;(AggressiveOrders3(1+n:end,6)-AggressiveOrders3(1:end-n,6))./AggressiveOrders3(1:end-n,6)]; %#ok<AGROW>
        BetasMech{filecounter,1} = betasMech;
        %BetasMech05d{filecounter,1} = betasMech05d;
        StandardErrorsMech{filecounter,1} = standardErrorsMech;
        %StandardErrorsMech05d{filecounter,1} = standardErrorsMech05d;
        BetasNMech{filecounter,1} = betasNMech;
        StandardErrorsNMech{filecounter,1} = standardErrorsNMech;


        LambdasNMechL{filecounter,1} = lambdasNMech;
        BetasNMechL{filecounter,1} = betasNMechL;
        SandwichesNMechL{filecounter,1} = sandwichesNMech;
        DaySizesNMechL{filecounter,1} = daysizesNMech;
        R2sNMechL{filecounter,1} = r2sNMech;


        LambdasMechL{filecounter,1} = lambdasMechL;
        BetasMechL{filecounter,1} = betasMechL;
        SandwichesMechL{filecounter,1} = sandwichesMechL;
        DaySizesMechL{filecounter,1} = daysizesMechL;
        R2sMechL{filecounter,1} = r2sMechL;
        VolumesMechL{filecounter,1} = volumesMechL;
        AdditionalValues{filecounter,1} = additionalValues;
        
        DatesForLater{filecounter,1} = datesForLater;
    end

    %Having finished calculating values for all 5 contract months,
    %Cycle through and store values in vectors with daily estimates
    betasExtrasMech = zeros(300,3);
    stErExtrasMech = zeros(300,3);
    counter = 0;
    pointer = 1;
    while(counter<5)
    counter = counter +1;
    lengthBetas = size(BetasMech{counter,1},1);
    betasExtrasMech(pointer:pointer+lengthBetas-1,:) = BetasMech{counter,1};
    stErExtrasMech(pointer:pointer+lengthBetas-1,:) = StandardErrorsMech{counter,1};
    pointer = pointer + lengthBetas;
    end
    betasExtrasMech = betasExtrasMech(1:pointer-1,:);
    stErExtrasMech = stErExtrasMech(1:pointer-1,:);
    
    betasExtrasMech05d = zeros(300,2*size(betaValues,1));
    stErExtrasMech05d = zeros(300,2*size(betaValues,1));
    
    betacounter = 0;
    while(betacounter<size(betaValues,1))
        betacounter = betacounter+1;
        counter = 0;
        pointer = 1;
        while(counter<5)
        counter = counter +1;

        lengthBetas = size(BetasMech05d{counter,1},1);
        betasExtrasMech05d(pointer:pointer+lengthBetas-1,2*betacounter-1:2*betacounter) = BetasMech05d{counter,betacounter};
        stErExtrasMech05d(pointer:pointer+lengthBetas-1,2*betacounter-1:2*betacounter) = StandardErrorsMech05d{counter,betacounter};
        pointer = pointer + lengthBetas;
        end
    end
    
    betasExtrasMech05d = betasExtrasMech05d(1:pointer-1,:);
    stErExtrasMech05d = stErExtrasMech05d(1:pointer-1,:);


    betasExtrasNMech = zeros(300,3);
    stErExtrasNMech = zeros(300,3);
    counter = 0;
    pointer = 1;
    while(counter<5)
    counter = counter +1;
    lengthBetas = size(BetasNMech{counter,1},1);
    betasExtrasNMech(pointer:pointer+lengthBetas-1,:) = BetasNMech{counter,1};
    stErExtrasNMech(pointer:pointer+lengthBetas-1,:) = StandardErrorsNMech{counter,1};
    pointer = pointer + lengthBetas;
    end
    betasExtrasNMech = betasExtrasNMech(1:pointer-1,:);
    stErExtrasNMech = stErExtrasNMech(1:pointer-1,:);

    lambdaExtrasNoMech = zeros(300,1);
    betaExtrasNoMech = zeros(300,2);
    sandwichesExtrasNoMech = zeros(300,4);
    daySizesNoMech = zeros(300,1);
    r2sNoMech = zeros(300,1);

    counter = 0;
    pointer = 1;
    while(counter<5)
    counter = counter +1;
    lengthLambda = size(LambdasNMechL{counter,1},1);
    lambdaExtrasNoMech(pointer:pointer+lengthLambda-1,1) = LambdasNMechL{counter,1};
    betaExtrasNoMech(pointer:pointer+lengthLambda-1,:) = BetasNMechL{counter,1};
    sandwichesExtrasNoMech(pointer:pointer+lengthLambda-1,:) = SandwichesNMechL{counter,1};
    daySizesNoMech(pointer:pointer+lengthLambda-1,:) = DaySizesNMechL{counter,1};
    r2sNoMech(pointer:pointer+lengthLambda-1,:) = R2sNMechL{counter,1};
    pointer = pointer + lengthLambda;
    end
    lambdaExtrasNoMech = lambdaExtrasNoMech(1:pointer-1,:);
    betaExtrasNoMech = betaExtrasNoMech(1:pointer-1,:);
    sandwichesExtrasNoMech = sandwichesExtrasNoMech(1:pointer-1,:);
    daySizesNoMech = daySizesNoMech(1:pointer-1,:);
    r2sNoMech = r2sNoMech(1:pointer-1,:);

    lambdaExtrasMech = zeros(300,1);
    betaExtrasMech = zeros(300,2);
    sandwichesExtrasMech = zeros(300,4);
    daySizesMech = zeros(300,1);
    r2sMech = zeros(300,1);
    volumesMech = zeros(300,1);

    counter = 0;
    pointer = 1;
    while(counter<5)
        counter = counter +1;
        lengthLambda = size(LambdasMechL{counter,1},1);
        lambdaExtrasMech(pointer:pointer+lengthLambda-1,1) = LambdasMechL{counter,1};
        betaExtrasMech(pointer:pointer+lengthLambda-1,:) = BetasMechL{counter,1};
        sandwichesExtrasMech(pointer:pointer+lengthLambda-1,:) = SandwichesMechL{counter,1};
        daySizesMech(pointer:pointer+lengthLambda-1,:) = DaySizesMechL{counter,1};
        r2sMech(pointer:pointer+lengthLambda-1,:) = R2sMechL{counter,1};
        volumesMech(pointer:pointer+lengthLambda-1,:) = VolumesMechL{counter,1};
        pointer = pointer + lengthLambda;
    end

    lambdaExtrasMech = lambdaExtrasMech(1:pointer-1,:);
    betaExtrasMech = betaExtrasMech(1:pointer-1,:);
    sandwichesExtrasMech = sandwichesExtrasMech(1:pointer-1,:);
    daySizesMech = daySizesMech(1:pointer-1,:);
    r2sMech = r2sMech(1:pointer-1,:);
    volumesMech = volumesMech(1:pointer-1,:);

    clearvars a2 AggressiveOrders3 dataIsNaNFormat diffjumps ESU1820180801 RegMtrx2 selector t2H t2M ;
    %plot(betasExtrasMech(:,3))

    toc


holidays = [11 28; 1 20; 2 17; 5 31; 7 3; 9 7];
selector = volumesMech>=100000;

%Store values
datesExtras = zeros(300,1);
counter = 0;
pointer = 1;
while(counter<5)
counter = counter +1;
lengthdates = size(DatesForLater{counter,1},1);
datesExtras(pointer:pointer+lengthdates-1,:) = DatesForLater{counter,1};
pointer = pointer + lengthdates;
end
datesExtras = datesExtras(1:pointer-1,:);

%Filter out holidays
selectorNonHolidays = zeros(315,1)+1;
selectorNonHolidays([8 63 127 192 257],:) = 0;
selectorNonHolidays = (selectorNonHolidays>0);
formatOut = 'mm/dd/yy';
dateVec = datestr(datesExtras+1,formatOut);

str = strcat('output\ComprehensiveStudy_6',Strcatroots{1,outerloop},'AllVarying_030through060');

%Save to file
save(str);

end

